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Cubic spline interpolation is a mathematical procedure which is 
an analog of the draftsman’s plastic spline. The advantage ofthis | 
interpolation procedure over the more commonly used methods 
such as Lagrange lies in its ability to not only fit each given data 
value exactly but to maintain continuity of the first and second de- 
rivatives. 


The relative accuracy of the cubic spline interpolation procedure 
for generating gridded data values and estimates of mean gravity 
anomalies from track-type geophysical surveys is shown to be ex- 
cellent when applied to properly designed surveys. Techniques for 
interpreting the two-dimensional Fourier transform in terms of track 
spacing, track orientation, and down-track sampling rate are pre- 
sented to demonstrate the effect of these parameters on interpolation 
accuracy. A procedure for utilizing closed form integration of the 
bicubic spline surface to produce mean gravity anomalies is shown 
to yield accuracies comparable to the method of averaging cubic 
spline grid values. 
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FOREWORD 


The U. S. Naval Oceanographic Office, as part of its mission to 
provide deflection of the vertical information in ocean areas, is currently 
engaged in studies concerning the accuracy of computing deflections of the 
vertical from gravity survey data. The deflection of the vertical calculation 
is the last step in a process involving many stages; consequently, the overall 
accuracy will be dependent on the accuracy of each individual stage. This 
report is concerned with the interim step of determining mean anomalies for 
use in the Vening Meinesz formulation. Utilization of the numerical method 
presented in this report is expected to result in a significant improvement in 
the accuracy of interpolating mean anomalies from original survey data. 
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INTRODUCTION 


Many of the analytical techniques which are currently in use in geophysics 
involve some form of mathematical operation on digitized survey data. The two- 
dimensional character of some of these techniques requires that the survey data, 
usually collected along tracks, be reduced to values of the surveyed parameters 
on an equally spaced grid. Examples of these techniques include upward and 
downward continuation, derivative calculations, regional trend removal, two- 
dimensional trend filtering, model studies, deflection of the vertical and geoidal 
undulation calculations, and automated contouring techniques. The interpolation 
procedures which are discussed in this report are applicable to any measurements 
obtained by track-type surveys. Emphasis, however, is placed on the problem of 
gridding gravity survey data and the determination of mean values of the gravity 


field within specified areas. 


In general there are two approaches available for gridding track-type 
survey data. One approach is that of fitting a least-squares surface, usually of 
polynomial form, to all of the available survey data within a predetermined area 
and using the resulting formula to generate the required grid values. Examples 
of this procedure, applied to irregularly spaced data, may be found in Crain and 
Bhattacharyya (1967) and Krumbein (1959). The advantages of this least-squares 
method are that it tends to operate as a low-pass filter, which reduces the effects 
of noisy or erroneous data and it can be applied directly to unequally-spaced 


data points. The chief disadvantages are: 


1. The wavelengths of the features which will be accurately fit is a 
function not only of the degree of the polynomial but also of the size of the 
area which is selected. This causes some difficulty in controlling the technique 


in a production type operation. 


2. The method is quite sensitive to round-off errors and, except in 
cases where orthogonalization techniques are used, can become unstable for 
high degree surfaces. Recent work by Wampler (1969) discusses the problems 


associated with inverting ill-conditioned matrices using standard procedures. 


3. The technique requires a rather large amount of computer time, 
although this problem can be minimized through the utilization of recursive 
linear-regression techniques (Fagin, 1964) when additional data points are 


added to the set. 


The second general approach which is used to grid survey data is based 
on interpolation techniques which fit each data value exactly. In this case, 
it is assumed that any large errors or noise in the data have been removed prior 
to application of the gridding process, and that the survey coverage is sufficient 
to minimize aliasing in the data. Historically, this approach has consisted of 
connecting each pair of data points in such a way as to partition the survey 
area into non-overlapping triangular sections (Rankin, 1963) and fitting a plane 
to each of these sections, or by partitioning the survey area into a preliminary 
rectangular grid and interpolating within each of these areas by an exact-fitting 
two-dimensional polynomial (Crain and Bhattacharyya, 1967). The principal 
disadvantage of these methods is that no conditions are imposed on the derivative 
of the interpolating functions. Except in very smooth areas, this results in the 
generation of non-realistic gradients. The problem is somewhat reduced by the 
use of a two-dimensional weighting function of the type developed by Shepard 
(1968). Shepard's technique generates the data value at a desired grid location 
by applying a weighting function based upon the distance and direction to the 
neighboring data points. In addition, estimates of the derivatives at each data 
point are included. This removes the requirement that the partial derivatives be 
zero at each data point thus allowing the interpolating function to have extrema 


anywhere and not just at the original data points. 


In order to overcome some of the difficulties associated with these 
existing techniques and to effectively utilize the computational efficiency 
afforded by track-type survey data, an interpolation procedure was developed 
by Bhattacharyya (1969) which makes use of cubic and bicubic splines. The 
advantage of the cubic spline method, which is a mathematical analog of the 
draftsman's plastic spline, lies in its ability to not only fit the observed data 
values but to maintain continuity of the first and second derivatives. This 
paper concerns the development and testing of a cubic spline technique as 
applied to the computation of mean gravity anomalies. By utilizing simulated 
survey data obtained from a model gravity field, it is shown that the spline 
technique will produce accurate estimates of the mean gravity anomalies 
provided the original survey is properly designed. The use of a two-dimensional 
Fourier transform for a solution to the survey design problem is also presented. 
DEVELOPMENT OF THE CUBIC SPLINE ALGORITHM FOR GRIDDING 
TRACK TYPE SURVEY DATA 

The algorithm presented here is based upon the technique developed 
by Bhattacharyya (1969) with some modifications to the spline formulation in 
order to use the boundary conditions and development presented by Pennington 


(1965). 
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are the values of the second derivative of f(x) at the given data points then 
for the interval x, Sx S% 41 We have 
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This is simply the equation for a straight line passing through the points 
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Integrating (1) once yields 
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where Cy is the constant of integration. 
Integrating (1) a second time yields 
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With the requirement that f(x) match the observed data at the points 
and , the constants c, and c, may be evaluated from equation (3). 
aimee 1 2 : 
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Solving equations (4) and (5) for the values of Cy and Cor and substituting these 


into equation (3), yields the following formula for the cubic spline in the 


interval x <x SX ay 
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eee OO L.(6) 
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In order to solve equation (6), the second derivatives (2) 21.44) 
must be determined at the end points of the interval. This is accomplished 


by requiring that the first derivative, equation (2), be continuous at each 


given data point, i.e., 


lim f(x) = f(x,). 


This means that when equation (2) is written for the intervals 
(x 447%) and (rays the two values obtained for f (x,) must agree. 
Thus, from equation (2), we have for the interval (ye 


f (x)= Ss ee aay + aa Nee Xi GGRINTS AG OE, (8) 


By equating (7) and (8), the formula for determining the values of the unknown 


second derivatives becomes 
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Now as k takes on the values 2—- m-1, equation (9) will produce m-2 equations 
in the m unknowns (z, Zp eee ze In order to solve this system, two additional 
equations are needed. These additional equations are obtained through 
utilization of an appropriate set of boundary conditions. There are several 
reasonable boundary conditions which may be selected. Bhattacharyya (1969) 
used the assumption that the second derivative is zero at the points (x, x) 

In the present development we will use the boundary conditions proposed by 
Pennington (1965). These conditions state that the second derivatives at the 
two end points are a linear extrapolation of those associated with the adjoining 


interval. Thus, from equation (1) we have 
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For k = 2 this reduces to 
tz, (+ +1) eT (11) 
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and for k =m-1 we have 
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Combining equations (11) and (12) with the m-2 equations obtained from (9) for 
k =2, 3...m-1, yields a set of m equations in m unknowns (24125 oe Zz). 
With the exception of one additional element in the first and last row, this set 
of equations forms a tridiagonal matrix which is easily solved using the standard 
Gauss reduction method. In matrix form, the set is given by AZ =B where 


A is mxm, Z is mxl and B is mxl thus, 
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Utilizing this solution in equation (6) yields an interpolation formula 
which fits each given data value exactly, has continuous first and second 
derivatives and is a simple cubic polynomial in x within the interval between 


each pair of data points. 


Since the cubic spline is a function of only one independent variable, 
the data obtained along a survey track must be adjusted to lie on a straight 
line. The interpolation may then be accomplished by utilizing distance along 


track as the independent variable. The procedure presented by Bhattacharyya 


(1969) is followed for this operation. The steps in this technique are as 


follows: 


1. Taking the data from one track at a time, the position of the data 
points are converted into x, y coordinates and a least-squares straight line is 
fitted to these locations. Since there is no statistical significance attached 
to this operation, either x or y may be considered the independent variable. 
The computer program listed in the appendix considers x, or equivalently 
longitude, as the independent variable. If the survey tracks happen to run 
exactly north-south, the program should be modified to consider y as the 
independent variable. With the coordinates of the data points given as 


(x. Yer i =1,N), the least squares line is f(x) = Ay +P ay x where 
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2. The perpendicular distance between the straight line and each data 
point is determined and utilized to project the points orthogonally onto the 
line with an adjusted data value determined by applying an estimate of the 
local gradient, assumed to be a function of distance only. If (ry) are the 


coordinates, and V, the value of the Kth data point, then the perpendicular 


k 


distance between this point and the least squares line is given by 


| Ae oy Ae SO) 
ee 


The mapped coordinates (erp) of this point on the least-squares straight line 


are given by 
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If B. is less than a predetermined pivot distance, the value associated with 


the data point is unchanged. If B. is greater than this pivot distance then 
the adjusted value Ve associated with the mapped coordinates (xe ryp) is 


computed by 
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with i = { a } as appropriate. 


In the computer program for the cubic spline algorithm contained in 
the appendix, the pivot distance is selectable via a control card. This 
distance is usually set equal to the maximum distance which one would be 
willing to move a data point without changing its value. In order to minimize 
the propagation of the error associated with the assumption that the gradient 
correction in equation (18) is independent of direction, continuous survey 
tracks which deviate appreciably from a straight line should be broken up 


into smaller segments with each segment treated as a separate track. 


The mapped coordinates and adjusted data values may now be considered 


as irregularly-spaced digital samples from a function with the independent 


variable being distance along track from some arbitrary starting point, and 


the dependent variable being of course the adjusted data values. 


3. Utilizing the mapped data, the cubic spline, equation (6), is 
determined for each track. This equation may then be used to interpolate 
data values at the intersections of the straight least-squares track lines and a 
set of parallel lines whose spacing is equal to the desired final grid spacing. 
If the direction of the survey tracks is predominately east-west then the 
direction of the set of parallel lines is north-south. Similarly, for north- 
south tracks, the lines are run east-west. The computer program, contained 
in the appendix, is designed to operate on tracks in any direction except 
exactly north-south. The direction of the set of parallel grid lines is controlled 
by the direction of track line number one. Since the track number designation 
is arbitrary, this feature allows the user to determine the desired orientation 
(N-S or E-W) of the parallel grid lines in order to obtain as many intersections 


as possible. 


4. The interpolated data values generated in step 3 may be regarded 
as unequally-spaced digital samples from a function whose independent variable 
is distance along each of these parallel lines. Application of the spline pro- 
cedure in this cross-track direction produces the final interpolated values at 
the desired grid points. If mean anomalies are desired, grid points are generated 
at one-half the final grid spacing and the resulting nine points are averaged to 
produce the mean value for each grid cell. 
THE BICUBIC SPLINE ALGORITHM FOR THE GENERATION OF 
MEAN ANOMALIES 

In the previous development (Step 4) the mean value of the function 
is obtained by averaging the nine interpolated values located within the region 
in which the mean is desired. The accuracy of this method will, of course, 


depend on the spacing of the points used in the averaging process relative to 


the frequency content of the measured field. We now consider an alternative 
approach of determining mean anomalies which utilizes the bicubic spline 
method of de Boor (1962). In this approach, a two-dimensional cubic poly- 
nomial surface is obtained by use of the bicubic spline, and the mean value is 
generated directly by integration. The algorithm presented here is based on 

de Boor's development, modified to utilize the preceeding cubic spline formula- 
tion and to compute mean anomalies by integration of the resulting bicubic 


surface. The development of the algorithm proceeds as follows. 


Assume that we have available a set of data points with coordinates 
(x.ry.), i=0...1, | =0...J, and with data value U.. at the intersections of 
| 


a rectangular grid covering an area R (Figure 1). 


FIGURE 1. THE Re GRID OF DATA VALUES FOR INPUT INTO THE 
BICUBIC SPLINE ALGORITHM. 


Utilizing this data set with an appropriate set of boundary conditions specified 


around the perimeter of R, a two-dimensional cubic polynomial of the form 


1] 


m 
U(x,y) 07) = Dye en Dag) CEs ae sb ane) 
m=0 n=0 
where x 12% S%; and We ESESY i 


may be determined for each R. . rectangle within R. In analogy with the cubic 

spline, the unknowns @ ule for this two-dimensional bicubic surface are le 

by requiring that the surface fit each data point exactly, i.e., R(x. "y )= 
U(x, y) dU(x,y) 

and that the partial derivatives seul and aye be es, In 


order to accomplish this, the first step is to compute estimates of the partial 


derivatives at each data point by utilizing the cubic spline formulation. Let 


au.. aU. a°U.. a, 
i lites Caciiem ay p i sii ~ Oxdy dy 


In order to compute Pi at each grid point, equations (9,11,12) are utilized 

to obtain values of the second derivative Z. from the U. data values along the 
nth row. Equation (8) is then used with these Z. to produce values for Pa 
i=0... 1-1 with Fe being computed from equation (7). This process is 
repeated for each row of input data to generate the full P.. matrix. In 
identical fashion, qi and af are computed by operating on columns of U.. 
and P. Respectively. Utilizing the given data values and these computed 


Dentin at each grid point, the unknown, @ HL , may be determined for 


each R rectangle covering R by solving the matrix equation 
Alc, 41) KA y.-y 4) = (ail) (20) 
TO Soe si ie wil mn 


The proof for this equation is given in de Boor (1962). 


Normally, the U.. are given on a square grid, in which case we let 
h=(x. - x. _1) =, - Y._y). The elements of the matrices in equation (20) 


are then 
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Substituting the @ a elements into equation (19) produces a two-dimensional 
cubic polynomial formula representing the data within each of the R.. subareas 


covering R. The mean value of the field (C), within the Re rectangle, is 


computed from equation (19) by 


ij 
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xX. 
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Si ne f 
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The bicubic spline computer program given in the appendix contains 
the option of evaluating this integral between arbitrary limits within each 
R.. rectangle. This feature allows computation of the mean anomalies at a 
grid spacing either equal to that of the U.. input data, or one-half or one-third 


the spacing of the input data. Setting the (x.y Y;-) origin equal to (0,0), 
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and letting the limits of integration be C< x <D and A<y <B, and with 


— = e e . 2 . 
co Lhe de h, integration of equation (21) yields 
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The computer program for the bicubic spline algorithm uses this formula to 
determine estimates of the mean anomaly on a specified grid interval. 
GENERATION OF MODELS FOR TESTING THE ACCURACY OF 
INTERPOLATION PROCEDURES 

Two essential criteria for determining the accuracy of an interpolation 
technique are (1) the existence of some reference surface with which the 
interpolated values can be compared, and (2) all errors due to sources 
unrelated to the interpolation process be minimal. We are concemed with 


testing interpolation methods for use with data collected by systematic surveys. 


For this specific application, the above criteria can be met by properly designed 


simulated surveys of a mathematically-derived reference surface. Proper design 


implies that the survey parameters viz., track spacing, track orientation, and 


(22) 


down-track sampling interval be chosen so as to define the essential characteristics 


of the reference surface. In the following material we discuss the theory and 


application of the fast Fourier transform for the design of surveys. In addition, 


the method of determining the reference surface for use in obtaining quantitative 


estimates of interpolation accuracy is given. 


A. The Two-Dimensional Fast Fourier Transform as a Tool for 
Simulated Survey Design 


Simulated surveys must be designed in such a way as to minimize 
aliasing of the spectral components of the data. At the present time, the 


sampling theory that we have available is based primarily on the Shannon 


14 


sampling theorem. This theorem states that a continuous band-limited function 
of time may be reconstructed from equally spaced digitized samples provided 
that you have at least two samples per shortest period. The reconstruction is 


accomplished by convolving the digital values with a Sinc function 


(2) ben 


Sin 20A (t - => 
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where the band limits are +2nA and the sampling rate is 2A. See Papoulis (1962) 
for a proof of this theorem. Although it is theoretically impossible for a 
function of finite length to be band limited, a reasonably accurate estimate 

of A may be obtained for a small area by application of a two-dimensional 


Fourier transform. 


In one dimension, the direct Fourier transform of an arbitrary function, 


g(x), is defined as 


(09) 
Gif) = af g(x) ou nes dx, (24) 
with f, equal to the frequency. By means of this formula, it is possible to 
transform information from the space domain into the frequency domain. G(f,) 
is, in general, a complex number of the form a(f, ) + ib(f,) or equivalently 
2 ely) 2ae 


b 
a +b We Sta a where ie ze |8) }) 


= |G) | is called the 
amplitude spectrum and arctan ae O(f,) is termed the phase spectrum. For 
our application of the Fourier transform, we will be primarily concerned with 
the amplitude spectrum. In the context in which it is used here, the term 


frequency refers to spatial frequency with the units being either cycles per unit 


distance or normalized to cycles per data interval . 


For a two-dimensional function g(x,y), equation (24) is expanded to 


Bradt p ~i2u(f x +f y) 
Gif aL) = f f a g(x by x" y” dxdy . (25) 
-@O -@ 


With g(x,y) consisting of digitized values on a uniform grid covering a finite 

area R, an estimate of G(F rf ) may be obtained for the normalized frequency 
range 0 < hank < 0.5 cycles/data interval, by use of the discrete two-dimensional 
fast Fourier transform (FFT) of Cooley and Tukey (1965). This discrete approxima- 
tion to Na must be interpreted to provide estimates of the three survey 
parameters; track spacing, sampling rate, and track orientation. The specifica- 
tion of track orientation requires a knowledge of what happens to trends or 
lineations in the (x,y) domain when they are mapped into the (Frf ) plane. 

It has been determined by Fuller (1967) that a trend in the (x,y) plane maps 

into a trend in the Co) plane in an orthogonal direction. Fuller shows this 


simply in the following way. 


Let g(x,0) be an arbitrary function of x defined for y =o. Using this 
definition, we can construct a function of both x and y, i.e., g(x,y), by 
shifting g(x,o) a distance Sy in the x direction, with S = ee as shown in 


Figure 2. 


FIGURE 2. THE FUNCTION g(x,y). 
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We will use what is usually called the "shift theorem" from transform theory 
which states that, if Gif) is the Fourier transform of g(x) then, from equation 
(24), eieal, x Gif) is the Fourier transform of g(x-A). With A = Sy, this 
theorem is used to integrate out the x dependence in equation (25) resulting in 


the equation 


inf SHE) 


Y 
iG) ll SGNe (26) 
= 

and, by symmetry, 

Ff) = 2G(F) fi os [2ny(F, s+] dy. (27) 
Integrating over y yields 

Sin [2nvir, S +) 
Ait) 280) saga 
x y 
FEA) = 2YG(E) Sine [2av (FS +f)] ; (28) 


Since the Sinc function has its maximum value when the argument is zero, 


equation (28) defines a function in the (fer f ) plane possessing a trend with 
-f 
Ss == . Thus we see that, since trends in a (x,y) plane are mapped into 
x 
trends in the (Ff ) plane in the orthogonal direction, survey track direction 


may be determined from the two-dimensional FFT by aligning the tracks parallel 

to the trends in the (fs e) plane. With this direction specified, the shape of 

the two-dimensional amplitude spectrum is used in conjunction with the Shannon 
sampling theorem to estimate the band limits of g(x,y) which, in turn, define the 
required track spacing and sampling rate for the survey. A general two-dimensional 


FFT computer program utilizing equally spaced gridded data and a fast Fourier 
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transform subroutine called NLOGN, written by Robinson (1967), is contained 
in the appendix. This program was used to design the simulated surveys for 


testing the cubic spline algorithms. 
B. Development of a Point-Mass Model Field 


To provide a reference surface for determining quantitative estimates of 
interpolation accuracy and to ensure that the data used in the interpolations are 
error free, gravity model data were generated from a random distribution of 


point masses as follows. 


Consider a point (P) on the surface of a sphere (Figure 3) with spherical 
coordinates (8, \), and a point mass located at depth d beneath the surface with 


coordinates (@', X'). The gravitational potential at P(®, X) due to the point mass is 


-km 
Vp =a (29) 


where m is mass, k is the universal gravitational constant, and r is the distance 


between the point mass and the point (P). 


P(®,)) 


FIGURE 3. GEOMETRY FOR COMPUTING THE GRAVITATIONAL 
POTENTIAL DUE TO A POINT MASS. 


For a sphere of radius R, the distance r is 


r= (Ro eR? - 2RR' cos Le 


where R'=R-d 
and 
cos W = cos® cos 9’ + Sin® Sin @' cos(X-dA). 


The radial (vertical component) Gp of the gravitational force at P is 


oN km 
Gp(9,4)= - R- = a (R-R' cos). 


For a finite number (N) of point masses situated at various depths and locations 
beneath the surface of the sphere, we can generate gravity anomalies on the 
sphere by summing the radial components due to the N point masses, i.e., 
Ne ne 
Gp (8A) = ZS gz (R-Ri cos p.). (30) 


i=l (F 
I 


A computer program to calculate the radial, easterly, and northerly components, 
in addition to the deflection of the vertical components with respect to a 
spherical earth of a set of point masses is given in the appendix. 
COMPARISON OF GRIDDING TECHNIQUES AND THE COMPUTATION 
OF MEAN ANOMALIES USING MODEL DATA 

The accuracy of the spline interpolation method, relative to several 
other techniques, is now considered in terms of gridding survey data and the 


computation of mean gravity anomalies by utilizing simulated survey data. 


A set of gravity anomalies, for use as a reference surface, were 
calculated from equation (30), at one-minute intervals, over an area of 80 x 


80 square minutes. The anomalies, given in Figure 4, result from 16 point 


| 


FIGURE 4. CONTOUR CHART OF THE GRAVITY ANOMALIES GENERATED 
BY 16 POINT MASSES WITH CONTOUR INTERVAL OF 5 mals. 


masses located at depths ranging from 7 to 12 kilometers beneath the surface 
of a spherical earth. Amplitudes of the anomalies range from -17 mgls to +83 
mgls. For purposes of comparing the various models and residual charts to 
follow, a transparency of Figure 4 is contained in a pocket attached to the 


rear cover of this report. 
A. Gridding of Survey Data 


To design a simulated survey for obtaining data to test the interpolation 
procedures, a two-dimensional FFT was computed using all of the model field 
data (Figure 4) available on a one-minute grid interval. Contours of the 
amplitude spectrum resulting from this computation are shown in Figure 5. 

Since only the shape of the spectrum is used to design the simulated survey, 

the contour values are arbitrarily normalized, and the frequency range is 
normalized to the units of cycles per data interval. Utilizing the result 
embodied in equation (28), Figure 5 shows that the trends in the model field 
are predominately in the NE-SW direction indicating that the preferred survey 
track direction is NW-SE. With the track direction specified, the flattening 
of the amplitude spectrum at a frequency of approximately 0.17 cycles per data 
interval in both the along-track and cross-track directions indicates that a track 
spacing and sampling interval of three nautical miles would be sufficient to 
minimize aliasing of the data. Figure 6 shows the simulated survey, designated 
as model 1, which was designed from this interpretation of the two-dimensional 
amplitude spectrum. By overlaying this survey with the transparent copy of 
Figure 4, it is apparent that no attempt was made to sample the peak values of 
the anomalies. Simulated deviation of the survey vehicle from the planned 
survey tracks was included, but measurement errors were excluded by computing 
the exact value of the gravity field at each survey point by using equation (30). 
The cubic spline algorithm was utilized with this survey data to compute inter- 


polated values of the gravity field at one minute grid intervals. The interpolated 
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FIGURE 5. THE TWO-DIMENSIONAL FOURIER TRANSFORM OF THE ENTIRE 
80 X 80 MINUTE MODEL FIELD WITH ONE DATA INTERVAL 
EQUAL TO ONE MINUTE AND CONTOURS IN NORMALIZED 
UNITS. 
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FIGURE 6. 


MODEL 1 - SIMULATED SURVEY TRACKS BASED ON THE 
FOURIER TRANSFORM OF THE 80 X 80 MINUTE MODEL 
FIELD. NOMINAL TRACK SPACING AND SAMPLING 
INTERVAL EQUAL TO 3 nm. 
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value at each grid point was then subtracted from the true value obtained via 
equation (30) and the result was machine contoured at a contour interval of 
0.5 mgl. This residual contour chart is shown in Figure 7. Figure 8 isa 


histogram of the model-1 errors with the standard deviation and range noted. 


In order to compare the accuracy of the spline interpolation to that 
which could be achieved through the use of a two-dimensional least-squares 
polynomial surface, a twelfth-degree surface was fitted to the simulated survey 
data obtained from model 1. The contoured field values generated by this 
surface are shown in Figure 9. Comparison of this surface with the true field 
values in Figure 4 indicates that the least-squares surface reproduces the 
general shape of the low-frequency components but does not depict the higher 
frequency features. As stated previously, the least-squares surface can be made 
to fit the high frequency components by reducing the size of the area over which 
the surface is computed but this is not considered to be an appropriate approach 


to interpolation when adequate survey data is available. 


By overlaying Figure 7 with the contour chart of the model field, it is 
readily apparent that the errors in the interpolated values are correlated with 
those regions of the model field which contain a large amount of energy at the 
higher frequencies. This correlation is interpreted as arising from the fact that 
the two-dimensional FFT is essentially a least-squares operation which produces 
an estimate of the average amplitude spectrum over the entire region in which 
the computation is carried out. This effect, coupled with the obvious non- 
stationarity (in space) of the model field, caused the amplitude spectrum to 
flatten out at an unrealistically low frequency (0.17 cycles/data interval). 

In turn, this led to an estimate of track spacing and sampling interval which was 
too large to adequately define the higher frequencies. In order to determine the 
effect of a shorter sampling interval, a second model, designated model 2, was 
constructed with track spacing and direction similar to model 1 but with a 


sample interval of 1 nm. The histogram of the errors resulting from the 
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FIGURE 7. CONTOUR CHART OF MODEL-1 POINT-VALUE RESIDUALS 
(SPLINE INTERPOLATION VS. TRUE VALUE) WITH A CONTOUR 
INTERVAL OF 0.5 mgl. 
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FIGURE 9. CONTOURS OF THE 12° LEAST-SQUARES POLYNOMIAL 


SURFACE GENERATED FROM MODEL-1 DATA. CONTOUR 
INTERVAL IS 5 mgls. (SURFACE COMPUTED FROM PROGRAM 
SUPPLIED BY GULF RESEARCH AND DEVELOPMENT COMPANY .) 
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interpolated values determined from this survey is shown in Figure 10. It is 
apparent that only a slight decrease in error was achieved by use of the increased 
sampling rate. In order to more effectively reduce this error, the FFT was 
computed for that anomalous portion of the model field contained within 10’ 

and 70' latitude and 0' and 60' longitude. Contours of the normalized amplitude 
spectrum for this part of the model field are shown in Figure 11. With an estimate 
of 0.25 cycles per data interval as the frequency at which the spectrum for this 
smaller area flattens out, a track spacing and sampling interval of 2 nm was 
selected. Although the lack of a definite trend in the high frequency components 
indicates that the track direction is not a particularly significant factor for this 
model field, the distinct trending of the low-frequency components indicates 

that some control of this factor is required. The model 3 simulated survey, based 
on the amplitude spectrum in Figure 11, is shown in Figure 12. Figure 13 isa 
contour of the differences between the spline interpolated data values obtained 
from the model-3 survey and the true values from the model field. The histogram 
of these errors, as shown in Figure 14, indicates that the model-3 survey produced 


a significant decrease in the overall errors in the interpolated data values. 


In order to determine whether the model-3 errors were caused by the 
cubic spline interpolation algorithm or by the simulated survey design, a two- 
dimensional, high-pass filter was applied to the one-minute grid values of the 
model field. This digital filter, which was designed by Fuller (1967) as a 
general high-pass filter, has the two-dimensional amplitude response shown in 
Figure 15. The weight function for this filter is symmetric, thus, no phase 
shifting is involved in the operation. The amplitude response is such that those 
frequencies which were not adequately sampled with the model-3 survey are 
retained while the lower frequencies are removed. The contoured result of 
this filtering operation is shown in Figure 16. By comparing this result with the 
residual contours in Figure 13, it is apparent that at least some of the errors in 
the interpolated values from the model-3 survey are generated by locally 


inadequate sampling. 
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FIGURE 11. THE TWO-DIMENSIONAL FOURIER TRANSFORM OF THE MODEL 
FIELD FROM LONGITUDE 0' TO 60' AND LATITUDE 10' TO 70' 
WITH ONE-DATA INTERVAL EQUAL TO ONE MINUTE AND 
CONTOURS IN NORMALIZED UNITS. 
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FIGURE 12, 


30’ 60’ 


SIMULATED SURVEY TRACKS BASED ON THE FOURIER TRANSFORM 
OF THE 60 X 60 MINUTE MODEL FIELD. NOMINAL TRACK 
SPACING AND SAMPLING INTERVAL EQUAL TO 2 nm. 
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FIGURE 13. CONTOUR OF MODEL-3 POINT-VALUE RESIDUALS (SPLINE 
INTERPOLATION VS. TRUE VALUE) WITH A CONTOUR 
INTERVAL OF 0.5 mgl. 
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FIGURE 15. AMPLITUDE RESPONSE OF THE TWO-DIMENSIONAL HIGH-PASS 
FILTER APPLIED TO THE 80 X 80 MINUTE MODEL FIELD. CONTOURS 
ARE IN NORMALIZED UNITS. 
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FIGURE 16. CONTOURS OF POINT-MASS ANOMALY DATA AFTER APPLICATION 
OF HIGH-PASS FILTER. CONTOUR INTERVAL IS 0.1 mgls. 
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These tests lead to the conclusion that, given an adequate track type 
survey design based on the frequency characteristics of the data, the cubic 
spline algorithm is a reliable interpolation technique. The computational 
efficiency of the algorithm is demonstrated by the fact that the 81 x 81 grid 
values interpolated from the model-1 survey required only one minute and ten 


seconds of computer time on a CDC-3800. 


As a qualitative test of the cubic spline gridding procedure applied to 
actual shipboard survey data, the algorithm was used to grid total magnetic 
intensity data obtained from a well-controlled marine magnetic survey consisting 
of E-W tracks spaced nominally 3 nm apart. Figure 17 is the hand-drawn 
contour chart of this area and Figure 18 is the machine contoured result using 
the cubic spline generated grid points. Aside from the fact that the cubic 
spline generated what is considered to be a more realistic position for the peaks 
of the anomalies, the two contour charts are nearly identical. The computer 
time required to generate the gridded data and plotter tape was less than three 
minutes on a CDC 3800 with an actual plotter time of twenty minutes. This is 
compared to an estimate of approximately 100 man-hours to plot the data values 


and produce the hand-drawn contour chart. 
B. Mean Gravity Anomaly Calculations 


In order to test the algorithms for computing mean gravity anomalies, a 
survey was designed which simulated a conventional shipboard gravity survey. 
This survey plan, designated model 4, consists of E-W tracks spaced 3 nm apart 
with a sampling interval of 1/3 nm. Figure 19 is a plot of this survey showing 
the random deviation from the planned track which was added to simulate a 
normal survey track. Figure 20 is a histogram of the errors in the interpolated 
data points resulting from the model-4 survey. Comparing this result with that 
obtained from the model-3 survey (Figure 14) indicates that the increased sampling 


rate (six times that of model 3) was sufficient to overcome the effect of increasing 
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FIGURE 19. MODEL 4 - SURVEY PLAN APPLIED TO POINT-MASS MODEL 
DATA TO SIMULATE A NORMAL SHIPBOARD GRAVITY SURVEY, 
i.e., TRACK DIRECTION E-W, TRACK SPACING 3 nm, 
SAMPLING INTERVAL 1/3 nm. 
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the track spacing and surveying at an angle to the trends in the data. Obviously 
this will not be true in the general case, since the effect of changing any of the 
survey design parameters will be controlled by the characteristics of the field 


being measured. 


Survey models 3 and 4 were used to produce data for testing three 
algorithms for computing mean gravity anomalies. The technique for generating 
mean anomalies using both a nine point averaging of the gridded data produced 
by the cubic spline and the two-dimensional integration of the bicubic spline 
surface have been discussed previously. In addition to these methods, a 
technique was tested which utilized a one-dimensional Lagrange interpolation 
(see e.g., Hamming, 1962). This procedure generated mean gravity anomalies 
on a uniform grid by averaging all of the survey data points falling within a 
grid interval to form a mean for that interval. Following this, a one-dimensional 
fourth-degree Lagrange interpolation was applied to these mean values to obtain an 
estimate of the mean gravity value for each of the grid cells for which no survey 
data was available. The results of these mean anomaly tests are shown in 
Figures 21-27. Figure 21 is a contour chart of the differences between the 
true mean anomalies computed for each one-minute square of the model field 
via equation (30), and the mean values computed by averaging the 1/2 minute 
grid values produced by applying the cubic spline procedure to the survey data 
from model 3. Figure 22 is a contour chart showing the differences between the 
true one-minute means and those generated by application of the Lagrange 
procedure to the model-3 data. The contour interval is the same in these two 
figures. The errors produced by the lack of control of the derivatives in the 
Lagrange interpolation method is apparent in Figure 22. It is probable that 
these errors could be reduced somewhat by using a third-degree Lagrange 
polynomial instead of the fourth-degree but, in light of the obvious advantages 


of the spline formulation, no additional testing of this approach was undertaken. 
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FIGURE 21. CONTOURS OF MEAN-ANOMALY RESIDUALS GENERATED 


FROM MODEL-3 SURVEY BY UTILIZATION OF A CUBIC 


SPLINE INTERPOLATION PROCEDURE WITH CONTOUR 
INTERVAL OF 0.5 mgl. 


FIGURE 22. CONTOURS OF MEAN ANOMALY RESIDUALS GENERATED 
FROM THE MODEL-3 SURVEY BY UTILIZATION OF A 


LAGRANGE INTERPOLATION PROCEDURE WITH CONTOUR 
INTERVAL OF 0.5 mg}. 
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Figure 23 compares the histograms of the errors generated by the cubic spline 
and Lagrange algorithms. Figures 24-26 show that simular results are obtained 
with the cubic spline and Lagrange methods when they were applied to a con- 
ventional shipboard gravity survey simulated by the test model 4 previously 


described. 


In order to determine if the increased computing time necessary for 
producing mean anomalies by integrating the bicubic spline surface was justified, 
the algorithm was applied to the one-minute grid values generated by application 
of the cubic spline to the model-3 survey data. A histogram of the errors 
resulting from this technique is shown in Figure 27. A comparison of Figures 
23 and 27 clearly indicates that the unique capabilities of the bicubic spline 


approach are wasted in so far as the computation of mean anomalies is concerned. 
CONCLUSIONS 


After evaluating the results obtained from this series of tests, the 


following conclusions may be drawn: 


1. The cubic spline algorithm for gridding track-type survey data and 
for computing mean gravity anomalies is accurate and computationaly efficient 


in comparison with the other techniques which were tested. 


2. The significant improvement achieved by the model-3 survey in 
comparison to the model-1 survey illustrates the concept of the sampling 
theorem embodied in equation (23), i.e., the accuracy of mathematical 
interpolation, and in particular the spline method, is highly dependent on 
the density of the survey data points relative to the frequency content of the 
anomalies. The results obtained by models 3 and 4 show that, if a survey is 
well designed, highly accurate interpolation may be anticipated by use of the 
spline procedure. Implicit in this conclusion is the requirement that the survey 


data itself be of high quality. 
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FIGURE 24, CONTOURS OF MEAN ANOMALY RESIDUALS GENERATED 
FROM MODEL-4 SURVEY BY UTILIZATION OF A CUBIC 
SPLINE INTERPOLATION PROCEDURE WITH CONTOUR 
INTERVAL OF 0.5 mgl. 
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FIGURE 25, CONTOURS OF MEAN ANOMALY RESIDUALS GENERATED 
FROM MODEL-~4 SURVEY BY UTILIZATION OF A LAGRANGE 
INTERPOLATION PROCEDURE WITH CONTOUR INTERVAL 
OF 0.5 mgl. 
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3. The results obtained by computing mean anomalies using a two- 
dimensional integration of the bicubic spline surface do not justify the 


additional computing time required for this technique. 
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APPENDIX A 


FORTRAN PROGRAM FOR GRIDDING AND COMPUTATION OF 
MEAN ANOMALIES USING THE CUBIC SPLINE 
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APPENDIX C 


FORTRAN PROGRAM FOR POINT-MASS MODEL DATA 
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